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Recent work has shown that the many-body expansion of the interaction energy can be used to develop 
analytical representations of global potential energy surfaces (PESs) for water. In this study, the role of 
short- and long-range interactions at different orders is investigated by analyzing water potentials that treat 
the leading terms of the many-body expansion through implicit (i.e., TTM3-F and TTM4-F PESs) and 
explicit (i.e., WHBB and MB-pol PESs) representations. It is found that explicit short-range representations 
of 2-body and 3-body interactions along with a physically correct incorporation of short- and long-range 
contributions are necessary for an accurate representation of the water interactions from the gas to the 
condensed phase. Similarly, a complete many-body representation of the dipole moment surface is found to 
be crucial to reproducing the correct intensities of the infrared spectrum of liquid water. 


I. INTRODUCTION 


It is known that the global potential energy sur¬ 
face (PES) of a system containing N interacting wa¬ 
ter molecules can be formally expressed in terms of the 
many-body expansion of the interaction energy as a sum 
over n-body terms with 1 < n < iV, 1 

V N {xi,.. .,x N ) = E V^ 1B \x a ) + J2v (2B) (x a ,x b ) 

a a>b 

+ E Vi3B \ x a,Xb,X c ) + . . . 

a>b>c 

+ V ( ' NB \x 1 ,...,x n ). (1) 

Here, Xi collectively denotes the coordinates of all atoms 
in the i-th water molecule, y( 1B ) is the one-body poten¬ 
tial describing the energy required to deform an individ¬ 
ual water molecule from its equilibrium geometry, and all 
higher n-body interactions, V^ nB \ are defined recursively 
through 

V^ nB \xx,... ,x n ) =V n (xi,...,x n ) - y^F (1B) (Zq) 

a 

-E^ ( 2 B| ( i «' i i’)- ■■■ 

a>6 

E y ((n - 1 ) B) ( a ; Qij;ra2; _ . . ,Xa n _ 1 ). (2) 

a,i>-">a n -i 

The rapid convergence of Equation (1), demonstrated in 
several studies of water clusters, 2-7 suggests that the PES 
associated with a system containing N water molecules 
can in principle be represented as a sum of low-order in¬ 
teractions that are amenable to accurate calculations us¬ 
ing correlated electronic structure methods [e.g., coupled 
cluster with single, double, and perturbative triple exci¬ 
tations, CCSD(T), which currently represents the gold 


standard in quantum chemistry]. 

The theoretical modeling of many-body effects in wa¬ 
ter began in the late 1960s when self-consistent field 
(SCF) calculations were carried out on small water 
clusters. 1,8-13 These studies concluded that nonadditive 
effects in water are generally nonnegligible, with 3B con¬ 
tributions being as large as 10 15% of the total pair 

interaction for ring structures. 4B effects were estimated 
to be, on average, ~1% of the total pair interaction. 
The first attempt to derive potential energy functions 
for water from ab initio data was made by Clementi and 
co-workers who fitted dimer energies calculated at the 
Hartree-Fock (HF) level of theory to an analytical rep¬ 
resentation of the 2B interactions. 14-16 Later, configu¬ 
ration interaction (Cl) calculations were used to derive 
the pairwise additive MCY potential, 17 to which 3B and 
4B terms were subsequently added through a classical 
description of polarization effects. 18 The first water po¬ 
tential including explicit 2B and 3B terms, derived re¬ 
spectively from 4 th order Moller-Plesset (MP4) and HF 
calculations, along with a classical description of higher- 
body polarization interactions, was reported in Ref. 19. 
Following these pioneering efforts, the ASP potentials 
with rigid monomers were derived using intermolecular 
perturbation theory. 20,21 

In the 2000s, Xantheas and co-workers introduced 
the TTM (Thole-type model) water potentials, which 
for the first time made use of a highly accurate IB 
term derived from ab initio calculations. 22-25 The lat¬ 
est versions (TTM3-F 26 and TTM4-F 27 ) were shown to 
reproduce the properties of water clusters, liquid wa¬ 
ter, and ice reasonably well, although some inaccura¬ 
cies were identified in the calculations of vibrational 
spectra. 28-30 Around the same time, 2B and 3B poten¬ 
tials with rigid water monomers were derived by Sza- 
lewicz and co-workers from symmetry-adapted pertur¬ 
bation theory (SAPT). 31,32 These studies eventually led 
to the development of the rigid-monomer CC-pol poten- 
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tial, a 25-site model with explicit 2B and 3B terms fitted 
to CCSD(T)-corrected MP2 dimer energies and SAPT 
trimer energies, respectively, with liigher-body terms be¬ 
ing represented through classical polarization. 33 CC-pol 
was shown to accurately reproduce the vibration-rotation 
tunneling spectrum of the water dimer and to predict the 
structure of liquid water in reasonable agreement with 
the experimental data. Within the CC-pol scheme, a re¬ 
fined 2B term with explicit dependence on the monomer 
flexibility has also been developed. 34,35 

More recently, the full-dimensional WHBB 36 and MB- 
pol 3 '~ 39 many-body potentials with flexible monomers 
were developed. WHBB includes explicit 2B and 3B 
terms fitted to CCSD(T) and MP2 data, respectively, 
with long-range many-body effects being represented by 
the same Thole-type model used in the TTM3-F poten¬ 
tial. MB-pol was derived from fits to CCSD(T) energies 
calculated for both water dimers and trimers in the com¬ 
plete basis set (CBS) limit and includes many-body ef¬ 
fects within a modified version of the polarization model 
employed by the TTM4-F potential. Although WHBB 
has successfully been applied to static calculations of sev¬ 
eral water properties, to date MB-pol is the only many- 
body potential that has been consistently employed in 
quantum molecular dynamics (MD) simulations that ac¬ 
curately predicted the properties of water from the gas 
to the condensed phase. 

While the calculations with many-body potentials re¬ 
ported in the literature demonstrate that Equation (1) 
can be effectively used to develop accurate molecular- 
level representations of the water interactions, the role 
of short- and long-range contributions at different orders 
has not been fully investigated. Here, we seek to ad¬ 
dress this issue by analyzing water potentials that treat 
the leading terms of the many-body expansion through 
implicit (e.g., TTM3-F and TTM4-F) and explicit (e.g., 
WHBB and MB-pol) representations, with a specific fo¬ 
cus on how these terms are defined and incorporated as 
a function of the intermolecular separations. The article 
is organized as follows: The four potentials are briefly 
described in section 2, an analysis of the energetics as 
well as of the vibrational frequencies and intensities for 
water systems ranging from the gas-phase dimer to the 
liquid phase is reported in section 3, and the conclusions 
are given in section 4. 


II. METHODS 

The four potentials studied here were chosen because 
they treat the intramolecular distortions on equal foot¬ 
ing, i.e., through the IB PES developed by Partridge 
and Schwenke. 40 In this way, the ability of the models 
to predict water properties is a direct reflection of the 
treatment of the intermolecular interactions. In addition, 
WHBB and MB-pol effectively share the same induction 
scheme employed by TTM3-F and TTM4-F, respectively, 
which thus allows for a systematic investigation of the 


effects of explicit low-order terms of the many-body ex¬ 
pansion. 

The dominant terms of the many-body water inter¬ 
actions are the 2B and 3B terms, which can be con¬ 
veniently split into short- and long-range contributions, 
y(2B.3B) = yJ h 2 ®’ 3B ) + V)f n g’ 3B) . At the two-body level, 
long-range interactions are dominated by electrostatics 
and dispersion, while exchange-repulsion and charge- 
transfer become increasingly important at short-range. 
3B interactions in water, on the other hand, arise primar¬ 
ily from induced interactions at long-range and exchange- 
repulsion when the monomer electron densities overlap. 41 
Both TTM3-F and TTM4-F include permanent and in¬ 
duced electrostatics as well as dispersion and repulsion at 
the 2B level, while all higher order terms are represented 
through many-body induction. 26,27 

Although WHBB and MB-pol seek to exploit the range 
separation of the low-order water interactions, the two 
potentials are intrinsically different both in philosophy 
and by construction. In WHBB, the short-range 2B 
term is described by a 7 th -degree permutationally invari¬ 
ant polynomial that smoothly transitions for an oxygen- 
oxygen separation between 6.5 and 7.5 A into the long- 
range component described by the 2B TTM3-F poten¬ 
tial. The 3B term only includes a short-range component, 
which is represented by either a 5 th - (WHBB5) or a 6 th - 
degree (WHBB6) permutationally invariant polynomial 
that dies off as the largest oxygen-oxygen separation be¬ 
tween two molecules of the trimer approaches 8.0 A. All 
higher-order interactions (n > 4) are described by the 
polarization model employed in the TTM3-F potential. 
By construction, WHBB thus employs a strict separation 
at the 2B and 3B levels between short- and long-range in¬ 
teractions that are described in a completely independent 
way and are essentially disentangled from all higher-order 
contributions. It should be noted that a simplified ver¬ 
sion of WHBB including only IB, 2B, and 3B interactions 
with shorter cutoffs was also used, 42 which, however, can¬ 
not be implemented in a straightforward way in standard 
software for molecular simulations in periodic boundary 
conditions. For this reason, only calculations with the 
full WHBB5 implementation of Ref. 36 are reported in 
the following analysis. 

On the other hand, MB-pol can be viewed as a classical 
polarizable model supplemented by short-range 2B and 
3B terms that effectively represent quantum-mechanical 
interactions arising from the overlap of the monomer elec¬ 
tron densities. Specifically, at all separations, lA 2B ) in¬ 
cludes (damped) dispersion forces derived from ab initio 
computed asymptotic expansions of the dispersion energy 
along with electrostatic contributions due to the inter¬ 
actions between the molecular permanent and induced 
moments. At short-range, V^ 2B ^ is supplemented by 
a 4 th -degree permutationally invariant polynomial that 
smoothly switches to zero as the oxygen-oxygen separa¬ 
tion in the dimer approaches 6.5 A. Similarly, the MB- 
pol 3B term, R( 3B ), includes a 3B polarization term at 
all separations, which is supplemented by a short-range 
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4 th -degree permutationally invariant polynomial that ef¬ 
fectively corrects for the deficiencies of a purely classical 
representation of the 3B interactions in regions where 
the electron densities of the three monomers overlap. 
This short-range 3B potential is smoothly switched off 
once the oxygen-oxygen separation between any water 
molecule and the other two water molecules of a trimer 
reaches a value of 4.5 A. In MB-pol, all induced interac¬ 
tions are described through many-body polarization us¬ 
ing a slightly modified version of TTM4-F. 37 MB-pol thus 
contains many-body effects at all monomer separations 
as well as at all orders, in an explicit way up to the third 
order and in a mean-field fashion at all higher orders. It 
should be noted that alternative functions to the permu- 
tationally invariant polynomials (e.g., GAP 43 ) have been 
suggested and successfully employed in modeling short- 
range many-body interactions. 44 

III. RESULTS AND DISCUSSION 

The effects of the different representations of the 
many-body interactions employed by TTM3-F, TTM4-F, 
WHBB, and MB-pol are investigated through the anal¬ 
ysis of the energetics of water systems ranging from the 
gas-phase dimer to the liquid phase. Figure la shows 
a comparison of the dimer binding energies calculated 
for the global minimum configuration of each PES. Al¬ 
though all potentials predict binding energies within ^0.2 
kcal/mol of the CCSD(T)/CBS reference value (inde¬ 
pendently of whether the dimer geometry is optimized 
within each PES or the CCSD(T)/CBS dimer geome¬ 
try is used in the calculations), 45 the analysis of the 
many-body water interactions reported in Ref. 46 shows 
that the TTM3-F and TTM4-F interaction energies devi¬ 
ate significantly from the CCSD(T) reference data when 
dimer configurations different from the corresponding 
equilibrium geometries are considered. Specifically, root 
mean square deviations (RMSDs) of 0.73 kcal/mol and 
0.44 kcal/mol were obtained in Ref. 46 for TTM3-F and 
TTM4-F, respectively, from comparisons with ~1,400 
CCSD(T)/aug-cc-pVTZ 2B energies corrected for the ba¬ 
sis set superposition error. It was instead shown that 
both WHBB and MB-pol provide a highly accurate rep¬ 
resentation of the whole multidimensional dimer PES 
as demonstrated by RMSDs of 0.15 kcal/mol and 0.05 
kcal/mol, respectively, obtained from comparisons with 
more than 40,000 CCSD(T)/CBS 2B energies. 37 

Due to the different treatment of many-body effects, 
the relative accuracies of the TTM3-F, TTM4-F, WHBB, 
and MB-pol potentials are expected to differ as the num¬ 
ber of molecules in the system of interest increases. This 
is shown in Figure lb, where the relative energies of low- 
lying hexamer isomers calculated with the four potentials 
are compared with the corresponding ab initio values. 47 
The hexamer isomers hold a special place in the exper¬ 
imental and theoretical studies of water because they 
are the smallest clusters in which the lowest energy con¬ 


figurations correspond to fully three-dimensional struc¬ 
tures. For this reason, the hexamer serves as a proto¬ 
typical model for the lrydrogen-bond networks observed 
in condensed phases. While TTM3-F, WHBB, and MB- 
pol (qualitatively) reproduce the isomer energy ordering, 
large deviations are found between the TTM4-F predic¬ 
tions and the reference ab initio values. Similar results 
are also obtained for other small water clusters, includ¬ 
ing the tetramer and pentamer isomers, shown in Fig¬ 
ures SI and S2 of the supplemental material. 49 . Since 
it was shown that the Thole-type polarization scheme 
employed by the TTM4-F potential correctly captures 
higher-order water interactions, 38,46 the low accuracy 
shown by TTM4-F in describing the relative energies of 
the hexamer isomers is possibly due to intrinsic deficien¬ 
cies in the 2B term. In this context, the performance 
of the TTM3-F potential, which provides the best agree¬ 
ment with the reference ab initio values, is somewhat sur¬ 
prising given the large deviations already seen at both the 
2B and 3B levels, 46 and may suggest some fortuitous can¬ 
cellation of errors. Based on this, it is interesting to note 
that the isomer energies predicted by WHBB, which em¬ 
ploys the same electrostatic model as TTM3-F, deviate 
from the reference data by as much as 2.0 kcal/mol. The 
origin of these deviations may result from inaccuracies 
in the short-range 2B and 3B WHBB polynomials, in¬ 
compatibility between the effective many-body represen¬ 
tation encoded in the TTM3-F electrostatic model with 
the short-range WHBB polynomials, or both. On the 
other hand, MB-pol, which employs a slightly modified 
version of the TTM4-F electrostatic model, predicts iso¬ 
mer energies in good agreement with the reference data, 
suggesting that the MB-pol short-range 2B and 3B terms 
accurately correct for the deficiencies that are intrinsic to 
the purely classical representation of these interactions 
encoded in the TTM4-F electrostatic model. 

To provide more quantitative insights into the differ¬ 
ent performance of the four potentials in describing the 
relative energies of the hexamer isomers, a many-body 
decomposition of the corresponding interaction energies 
(Table I) was carried out at the CCSD(T)-F12b/VTZ- 
F12 level of theory 50,51 with the VTZ-F12 basis set, 52 
which was shown to effectively provide the same accu¬ 
racy as CCSD(T)/CBS calculations at a reduced compu¬ 
tational cost. 53,54 As has been well established, 2 ~' the in¬ 
teraction energies are dominated by the 2B and 3B terms, 
with 4B and higher terms making largest (but still rel¬ 
atively small) contributions for cyclic structures. 55 It is 
interesting to note that, while the hexamer binding en¬ 
ergies predicted by TTM3-F are in excellent agreement 
with the reference calculations (Figure lb), the results 
of Table I demonstrate that this agreement is clearly 
achieved by a fortuitous cancellation of errors. For ex¬ 
ample, in the prism isomer, the TTM3-F 2B interaction 
is too attractive by ~3.2 kcal/mol while the 3B inter¬ 
action is too repulsive by ^3.8 kcal/mol, resulting in a 
total interaction energy that differs from the reference 
value by merely 0.79 kcal/mol. In contrast, TTM4-F, 



4 



FIG. 1. a) Binding energies for the global minimum configuration of the water dimer (in kcal/mol) obtained with the TTM3-F, 
TTM4-F, WHBB, and MB-pol potentials in comparison to the reference CCSD(T) value from Ref. 45(fillod symbols). Also 
shown as open symbols are the binding energies calculated using the reference CCSD(T) dimer geometry, b) Comparison of 
the relative binding energies (AE) of the water hexamer isomers calculated with the TTM3-F, TTM4-F, WHBB, and MB-pol 
potentials using the ab initio geometries of Ref. 47. Also shown as a reference are the ab initio values from Ref. 47. c) Mean 
absolute difference in the total energy between QMC and DFT with various exchange-correlation functionals and the TTM3- 
F, TTM4-F, WHBB, and MB-pol potentials calculated for configurations (in periodic boundary conditions) extracted from 
path-integral molecular dynamics simulations of water performed with the vdW-DF and vdW-DF2 functionals. The statistical 
errors in the mean absolute differences arising from the QMC calculations are less than 0.006 mHa/molecule. See Ref. 48 for 
specific details on the QMC and DFT calculations. 


which accurately represents the electrostatic interactions 
(as shown by 3B contributions that he within 1 kcal/mol 
of the reference values for all low-lying hexamer isomers), 
is unable to correctly describe the 2B interactions, which 
have significant contributions from exchange/repulsion 
and charge transfer. Thus, without the benefit of for¬ 
tuitous cancellation of errors (that can, to some extent, 
be encoded into models through empirical parametriza- 
tion), potentials with accurate electrostatics but without 
a comparably accurate treatment of short-range quan¬ 
tum mechanical effects may be expected to have rela¬ 
tively poor performance. Having been explicitly derived 
from many-body expansions of the interaction energies 
obtained from correlated electronic structure data, both 
WHBB and MB-pol closely reproduce the CCSD(T)-F12 
reference many-body terms. However, nonnegligible dis¬ 
crepancies in the 3B and 4B terms result in WHBB being 
appreciably less accurate than MB-pol at reproducing the 
hexamer interaction energies. 

The differences between the four potentials become 
more evident when their accuracy is assessed against 
benchmark quantum Monte Carlo (QMC) interaction en¬ 
ergies calculated for liquid water configurations extracted 
from path-integral molecular dynamics simulations of the 
vdW-DF and vdW-DF2 functionals. 48 QMC has been 
shown to be a reliable benchmark in the study of small 
water clusters, producing relative energies with an accu¬ 
racy comparable to that of CCSD(T). 58,59 As a measure 
of accuracy, the mean absolute deviation (MAD) between 
the energies £) obtained with each of the four potentials 
and the reference QMC energies i?Q MC was calculated as 


i Nc 

MAD = t- I ( E i - ^ QMC ) -(E- E qmc ) 

-/V c . I \ / 


( 3 ) 


Here, N c is the number of water configurations and 
(E — E^ mc ) = (1 /N c ) ~ E® MC ) is the average 

energy difference for all configurations, effectively align¬ 
ing the zero of energy with the reference QMC data. Fig¬ 
ure lc shows that the mean absolute deviation associated 
with WHBB is ~4 times larger than the corresponding 
values obtained with both TTM3-F and TTM4-F, and 
~15 times larger than the MB-pol value. Interestingly, 
MB-pol also achieves better accuracy than all DFT mod¬ 
els analyzed in Ref. 48. The TTM3-F and TTM4-F re¬ 
sults indicate that effective representations of the many- 
body interactions through classical polarization models 
can provide a reasonable description of the energetics as¬ 
sociated with the liquid phase (i.e., comparable with that 
provided by most of the DFT models currently used in 
water simulations), albeit through empirical parameteri¬ 
zation. It should be noted that, as shown by the results 
reported in Table I and the many-body analysis of dif¬ 
ferent DFT models carried out in Ref. 60, cancellation 
of errors between different terms of the many-body ex¬ 
pansion of the interaction energies may also affect the 
energetics of bulk systems described by both DFT and 
classical polarizable models. On the other hand, the in¬ 
creased accuracy of MB-pol relative to TTM4-F indicates 
that short-range many-body interactions beyond a purely 
classical electrostatic representation are necessary for a 
correct, molecular-level description of the liquid phase. 
The significantly different accuracies associated with the 
WHBB/TTM3-F and MB-pol/TTM4-F potentials also 
suggest that the correct integration of explicit short- 
range and effective long-range many-body interactions is 
critical for ensuring the transferability of the potential 
across different phases. 

The ability of TTM3-F, TTM4-F, WHBB, and MB- 
pol to represent the global multidimensional PESs of wa- 
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TABLE I. Many-body decomposition of the interaction energy of low-lying hexamer isomers. All energies are in kcal/mol. The 
CCSD(T)-F12b 50,51 calculations were performed with the VTZ-F12 basis set 52 and corrected for basis set superposition error 
through the site-site function counterpoise method 56 using the Molpro quantum chemistry package. 5 ' 


(a) Prism 



CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-38.65 

-41.83 

-34.96 

-38.59 

-38.96 

3B 

-8.78 

-5.02 

-8.41 

-9.19 

-8.74 

4B 

-0.66 

-0.42 

-0.48 

-0.42 

-0.52 

•5B 

0.06 

0.03 

0.06 

0.03 

0.05 

6B 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

-48.03 

-47.24 

-43.79 

-48.17 

-48.17 



(c) Book-1 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-35.76 

-40.22 

-35.21 

-35.53 

-35.80 

3B 

-10.46 

-5.65 

-10.05 

-9.99 

-10.26 

4B 

-1.08 

-0.56 

-0.89 

-0.56 

-0.92 

5B 

-0.04 

-0.02 

-0.01 

-0.02 

-0.01 

6B 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

-47.34 

-46.46 

-46.16 

-46.10 

-47.00 



(e) Bag 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-35.01 

-39.71 

-34.99 

-34.89 

-35.24 

3B 

-10.43 

-5.42 

-9.60 

-9.86 

-10.15 

4B 

-1.16 

-0.59 

-0.86 

-0.59 

-0.90 

•5B 

-0.02 

-0.02 

0.00 

-0.02 

-0.01 

6B 

0.01 

0.00 

0.00 

0.00 

0.00 

Total 

-46.61 

-45.74 

-45.45 

-45.36 

-46.30 



(g) Cyclic 

-boat-1 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-32.02 

-37.54 

-32.33 

-31.79 

-32.02 

3B 

-11.44 

-6.22 

-10.90 

-10.35 

-11.30 

4B 

-1.63 

-0.89 

-1.33 

-0.89 

-1.35 

•5B 

-0.16 

-0.09 

-0.09 

-0.09 

-0.09 

6B 

-0.01 

-0.01 

0.00 

-0.01 

0.00 

Total 

-45.26 

-44.75 

-44.65 

-43.13 

-44.76 


(b) Cage 


CCSD(T)-F12 TTM3-F TTM4-F WHBB5 MB-pol 


2B 

-38.20 

-42.01 

-36.43 

-37.99 

-38.48 

3B 

-9.04 

-4.77 

-8.92 

-9.45 

-8.93 

4B 

-0.53 

-0.26 

-0.43 

-0.26 

-0.47 

5B 

0.01 

0.01 

0.03 

0.01 

0.03 

6B 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

-47.77 

-47.02 

-45.74 

-47.69 

-47.85 



(d) Book-2 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-35.86 

-40.35 

-35.61 

-35.65 

-35.93 

3B 

-10.18 

-5.28 

-9.75 

-9.71 

-10.03 

4B 

-1.00 

-0.49 

-0.81 

-0.49 

-0.85 

5B 

-0.02 

-0.01 

0.00 

-0.01 

0.00 

6B 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

-47.06 

-46.14 

-46.17 

-45.86 

-46.81 



(f) Cyclic-chair 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-32.42 

-37.85 

-32.43 

-32.06 

-32.44 

3B 

-11.87 

-6.66 

-11.13 

-10.71 

-11.64 

4B 

-1.78 

-1.03 

-1.43 

-1.03 

-1.45 

5B 

-0.19 

-0.12 

-0.10 

-0.12 

-0.10 

6B 

-0.01 

-0.01 

0.00 

-0.01 

0.00 

Total 

-46.27 

-45.66 

-45.09 

-43.92 

-45.63 



(h) Cyclic 

:-boat-2 




CCSD(T)-F12 TTM3-F 

TTM4-F 

WHBB5 

MB-pol 

2B 

-31.96 

-37.55 

-32.45 

-31.66 

-32.02 

3B 

-11.43 

-6.35 

-10.93 

-10.45 

-11.29 

4B 

-1.61 

-0.90 

-1.33 

-0.90 

-1.35 

5B 

-0.16 

-0.09 

-0.09 

-0.09 

-0.09 

6B 

-0.01 

-0.01 

0.00 

-0.01 

0.00 

Total 

-45.16 

-44.91 

-44.81 

-43.11 

-44.76 


ter systems with an increasing number of molecules is 
directly reflected in the accuracy with which the four 
potentials predict the associated vibrational frequencies. 
Harmonic frequencies for the water hexamer have been 
reported at the CCSD(T)/aug-cc-pVDZ level. 62 More 
recently, ab initio reference harmonic frequencies for 
small water clusters have been obtained through two- 
body: many-body CCSD(T):MP2 calculations using the 
cc-pVQZ and aug-cc-pVQZ basis sets for H and O atoms, 
respectively, which enables a more complete treatment 
of the basis set through the hybrid two-body:many-body 
approach. 61 As shown in Figure 2, the WHBB and MB- 
pol PESs result in harmonic frequencies for the water 
dimer that deviate by less than 10 cm -1 from the ab ini¬ 
tio frequencies reported in Ref. 61. These results are con¬ 
sistent with the analysis of the dimer vibration-rotation 
tunneling spectra of WHBB and MB-pol, both of which 
exhibit nearly perfect agreement with the corresponding 
experimental data. 37 In contrast, relatively large devi¬ 
ations are seen for vibrational frequencies of both the 
TTM3-F and TTM4-F potentials. The different perfor¬ 
mance of the four potentials in predicting the relative 


energies of the hexamer isomers clearly leads to different 
levels of agreement with the ab initio harmonic frequen¬ 
cies. Independently of the isomer, the MB-pol values 
consistently lie within 50 cm -1 of the reference values 
while the WHBB harmonic frequencies can deviate, in 
some cases, by more than 100 cm -1 . Substantially larger 
deviations, up to ~200 cm -1 , are instead obtained with 
both TTM3-F and TTM4-F, reinforcing the notion that 
purely classical representations of the many-body water 
interactions are likely not sufficient to fully reproduce the 
complexity of the underlying multidimensional PESs. In¬ 
terestingly, all potentials predict somewhat larger devi¬ 
ations for the cyclic-chair isomer of the water hexamer, 
supporting early observations that many-body effects in 
water are relatively more important for ring structures. 13 

While the differences in the PESs clearly affect the 
underlying vibrational structure, as shown in Figure 2, 
the infrared activity of those vibrations is ultimately dic¬ 
tated by the associated dipole moment surfaces (DMSs). 
Many-body representations of higher-order electric prop¬ 
erties for molecular systems were characterized begin¬ 
ning in the early 1980s. 64 The first analysis of many- 
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Frequency (cm' 1 ) 

FIG. 3. Decomposition of the IR spectrum obtained from 
CMD trajectories with the MB-pol PES in terms of the many- 
body components of the MB-/r DMS. lB-Dip indicates that 
the one-body (gas-phase monomer) dipoles were used to cal¬ 
culate the dipole of the molecules sampled along the MB-pol 
CMD trajectories, from which the IR spectrum was calcu¬ 
lated. (1B+2B)-Dip indicates that short-ranged two-body 
dipoles were used in addition to the one-body dipoles. MB- 
Dip is the full MB -fi many-body dipole. The spectra were 
smoothed to facilitate the comparison between the line shapes 
obtained using the different approximations. The experimen¬ 
tal IR spectrum was taken from Ref. 63. 


body effects on the dipole moment of polar molecules 
was reported by Skwara et al. 65 As noted in 1996 by 
Schwenke, 66 molecular dipole moments can efficiently 
be represented in terms of geometry-dependent effective 
charges multiplied by their Cartesian positions. In line 


with these ideas, the first representation of the 2B dipole 
of water was reported as part of the WHBB suite. 36 A 
refined version of this IB + 2B DMS 42 has been used 
to calculate the transition dipole moments which were 
then used to modulate the WHBB5 frequency distribu¬ 
tions calculated within the local monomer approximation 
from configurations of liquid water extracted from MD 
simulations with the rigid E3B. 6 ' Neglecting both 3B 
and higher-order contributions to the dipole moments 
and dynamical effects (e.g., motional narrowing), good 
agreement in the relative intensities was obtained with 
the measured IR spectrum of liquid water. 42 Recently, a 
complete many-body representation of the DMS of wa¬ 
ter (MB-/z) was reported, 68 consisting of the one-body 
DMS of Lodi et al. 69 , an explicit two-body term, and 
a slightly modified version of TTM4-F polarization for 
long-range two-body and all higher-order many-body in¬ 
duced dipole moments. The particular functional form 
of MB-/i was derived from a systematic analysis of the 
many-body convergence of the electrostatic properties of 
water. 46 

To investigate the effects of many-body dipole mo¬ 
ments on the IR spectrum of liquid water, many-body 
molecular dynamics (MB-MD) simulations, 68 within the 
centroid molecular dynamics (CMD) formalism, were 
carried out with the MB-pol PES in combination with 
the MB-p, DMS. As shown in Figure 3, while the overall 
shape of the IR spectrum is captured in the IB + 2B 
representation of the DMS, 3B and higher-order contri¬ 
butions to the dipole moment significantly affect the IR 




























































































7 


TABLE II. Cost per molecular dynamics step for different PESs relative to q-TIP4P/F (a non-polarizable model with 3 
point charges per molecule and Lennard-Jones interactions between oxygen atoms). The system examined contains 256 water 
molecules in periodic boundary conditions. Timings are presented for TTM3-F and TTM4-F, the underlying electrostatics 
model used by WHBB and MB-pol, respectively. The additional cost of WHBB and MB-pol beyond their baseline electrostatics 
represents the computational cost of the short-range 2B/3B polynomials. WHBB is implemented as described in Ref. 70 using 
the WHBB polynomial libraries provided by the authors. All timings were performed in a modified version of DL_POLY2 using 
a single core of a typical Intel Xeon E5-2640 based workstation. 


Model 

Description 

Cost per MD step relative to q-TIP4P/F 

q-TIP4P/F 

Point charge 

l.Ox 

TTM3-F 

1 polarizable site/molecule 

7.3x 

TTM4-F 

3 polarizable sites/molecule 

8.5x 

WHBB 

TTM3-F electrostatics + 
empirical 2B dispersion + 

2B short-range 7 tft -degree polynomial + 
3B short-range 5 tft, -degree polynomial 

29,000x 

MB-pol 

TTM4-F electrostatics + 
ab initio 2B dispersion + 

2B short-range 4 tft, -degree polynomial + 
3B short-range 4 (,l -degree polynomial 

47x 


intensities. 2B contributions to the DMS lower the bend 
intensity while they contribute to ^70% of the intensity 
of the OH stretch band. Importantly, while 2B contribu¬ 
tions are critical to appearance of the shoulder at ~180 
cm -1 corresponding to the hydrogen-bonding stretch, the 
IB contributions are non-negligible. This indicates that 
some rotational motion is also involved in the hydrogen¬ 
bonding stretch. Nevertheless, the absolute intensities of 
the librational, bending, and stretching bands are only 
recovered when all many-body effects are included in the 
calculation of the dipole moment. 30 These results support 
early observations based on molecular dynamics simula¬ 
tions with the SPC potential supplemented with a 3B 
dipole-induced dipole term showing that the calculated 
far infrared spectrum of liquid water was in better agree¬ 
ment with experiment relative to results obtained includ¬ 
ing only a 2B description of the dipole moment. 71 

While using a highly accurate PES is a prerequisite 
for a physically correct description of the water proper¬ 
ties at the molecular level, an appropriate balance be¬ 
tween accuracy and efficiency is often critical when de¬ 
ciding which potential to employ in computer simulations 
since the computational cost directly affects the ability 
to calculate statistically converged quantities. As shown 
in Table II, a performance analysis carried out on a sin¬ 
gle Intel Xeon E5-2640 processor for a system consisting 
of 256 water molecules in periodic boundary conditions 
indicates that MB-pol is able to achieve a high level of 
accuracy at a cost of 47x that of the fixed point charge 
q-TIP4P/F model 72 and ~5.5x that of the TTM4-F po¬ 
tential. WHBB, on the other hand, is ~29,000 more 
expensive than q-TIP4P/F. A reference implementation 
of MB-pol is available as an independent plugin 73 for the 


OpenMM toolkit for molecular simulations. 74 


IV. CONCLUSION 

The role of short- and long-range contributions to the 
total energy of water systems ranging from the gas-phase 
dimer to the liquid phase was investigated by consider¬ 
ing four water potentials that treat the leading terms of 
the many-body expansion through implicit (i.e., TTM3-F 
and TTM4-F PESs) and explicit (i.e., WHBB and MB- 
pol PESs) representations. The analysis of the energet¬ 
ics, vibrational frequencies, and infrared intensity indi¬ 
cates that explicit short-range representations of 2B and 
3B interactions along with a physically correct incorpo¬ 
ration of short- and long-range contributions are neces¬ 
sary for an accurate representation of the water inter¬ 
actions, independently of the system size. These results 
thus suggest that atomistic water potentials built upon 
the many-body expansion of the interaction energy de¬ 
rived from ” first principles” hold great promise to achieve 
the long-sought-after goal of describing the macroscopic 
properties of water across different phases from a rig¬ 
orous microscopic viewpoint. A question that remains 
to be addressed is the extent to which these ’’first prin¬ 
ciples” many-body water potentials can be extended to 
the modeling of complex aqueous solutions. 
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